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Transition to turbulence in an elliptic vortex 


By T. S. Lundgren 1 AND N. N. Mansour 


1. Motivation and objectives 

We study the three dimensional instability and nonlinear growth of the two- 
dimensional flow described by the streamfunction 

_ Asmb l xsinb 2 y 
b\+b\ 

where b\ = n/Li, b 2 = 7r / L 2 . This is a swirling flow in a box which is bounded by 
0<x<Li,0<y<£2 and is infinite in the z direction. This flow is a solution 
of the Navier-Stokes equation with A = exp( — i/(6j + b\)t) which slowly decays. 
We seek a viscous solution which starts near this one and slips along but does not 
penetrate the bounding walls. The vorticity of the basic flow is ui z = A sin b\ x sin b 2 y 
which has maximum value A at the center of the box and drops to zero at the 
boundaries. We can think of the resulting flow as that of a captive vortex. 

Denote E = L\jL 2 . What is interesting about this flow is that when E ^ 1, 
it is unstable to three-dimensional disturbances. Periodic waves grow along the 
z direction, causing the captive vortex to distort into a snake-like configuration. 
Components of vorticity perpendicular to the z axis grow. As the instability grows, 
we find that the vorticity becomes increasingly sheet-like. Ultimately, the sheets 
break up into a turbulent mixture of intense vortex tubes. When E = 1, the flow 
is completely stable. 

The streamline pattern of the basic flow is elliptical in the central part of the box, 
with aspect ratio E , becoming more rectangular as the boundaries are approached. 
Figure 1 shows the streamline pattern when E = 2. Bayly (1989) has shown analyt- 
ically and numerically that the flow is unstable to high wavenumber disturbances 
when the Reynolds number is large. At high wavenumber, the instability in this 
bounded geometry is similar to the broad band instability in an unbounded flow 
with elliptical streamlines and uniform vorticity studied by Pierrehumbert (1986), 
Bayly (1986), Waleffe (1989, 1990), and Landman and Saffman (1987). The mecha- 
nism for the instability is a resonance between inertial oscillations, which can exist 
when the streamlines are circular, and the periodicity introduced by the ellipticity. 

The work which motivated these studies was that of Orszag and Patera (1983). 
They numerically found a three-dimensional instability on a flow which is a super- 
position of a Tollmein-Schlichting wave and the Blasius boundary layer flow. This 
flow contains an elongated elliptically shaped swirling region. It is believed that the 
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FIGURE 1. Streamline pattern and vorticity contours for E = 2. The tick marks 
are at .6 intervals in the r-direction and .4 intervals in the y-direction. 

wavy instability which develops explains the generation of streamwise vorticity in 
boundary layer transition. 

An even earlier work of Gledzer et al. (1974, 1975) showed instability in a swirling 
flow within an elliptic cylinder. The experiment consisted of rotating a water-filled 
elliptic cylinder until solid body rotation was achieved. Upon stopping the rotation, 
the water continues to rotate with approximately elliptical streamlines. It was 
observed that the flow rapidly developed swirls perpendicular to the rotation axis, 
with one or more cells depending on the length of the cylinder. Malkus (1989) did 
clever experiments with a water-filled flexible cylinder which was made elliptical by 
rotating it between stationary rollers. In this flow, the unstable wave rapidly flashes 
into small scale turbulence. 

There is another class of flows which is closely related to the present work. This 
is the bending instability of a sharp-edged vortex with uniform vorticity exposed to 
a transverse straining flow treated by Widnall, Bliss and Tsai (1974), Moore and 
Saffman (1975), and Robinson and SafFman (1984). In this flow, wavy disturbances, 
comparable in length to the dimensions of the vortex, grow on the vortex. In the 
absence of the straining flow, non-rotating planar waves exist. The strain causes 
such a wave to align with and grow along the plane of maximum strain rate. In our 
flow, the vorticity distribution rises smoothly to a maximum at the center instead 
of being sharp-edged, and there are finite boundaries; nevertheless, some of the 
stability results are very similar. 

In the present work, we show numerically that the basic flow described by Eq. 
(1) is unstable at modest axial wavenumbers when the Reynolds number is greater 
than a critical value. In section 2, we present stability diagrams which superficially 
resemble those for two-dimensional parallel flows. In section 3, we present the 
results of a nonlinear computation which clearly shows transition to turbulent flow. 

2. Stability 

For the stability problem, we set A = 1, neglecting the slow decay. This is similar 
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to the traditional method of studying the stability of a growing boundary layer by 
treating it as a parallel flow. In the usual way, we decompose the velocity field into 


the velocity field U of the basic flow plus a small perturbation u / . 
equations for u' are 

The linearized 

<i 

II 

o 

(2) 

du* 

— + V . (Uu' + u'U) = - Vp' + ^V 2 u'. 
ot 

(3) 

We separate variables in the form 


1 °° 

u ~expik z z — u(m^n)expi(rnb\x + nb 2 y) 

* m ,n=oc 

(4) 

1 °° 

v 9 =expik z z — \ 0(m, n)exp i(m6ix + nfe 2 y) 

bo ‘ ^ 

m,n=oc 

(5) 

oo 

w* = exp ik z z ^ u>(m,rc)expi(m&i:r + n^y) 

(6) 

m,n= oc 


with a similarly defined pressure variable, (Note that we have suppressed the time 
dependence in u, v and w.) Upon equating coefficients of exp i(k z z + mbix + nb^y) 
and using the continuity equation to eliminate the pressure, we get the following 

equations: 


du(m,n) bib?R x {m, n) 2l2 , 2t2 , l2x 

dt = b 2 +b 2 v \ rn Oj + n b 2 + k z ) u(m. 

.«) (7) 

dv(m, n) bifaRyim, n) , 2l2 , 2l2 _ l2x 

di = b 2 + b 2 ~ V ^ m +n b 2 + k z) «(m, 

n) (8) 

where 


d / \ rt \ rnb\(mF + nG + k z H) 

RAm,r,) = F(m,«)- ^ + ^ + ^ 

(9) 

p x n&^(mF + nG + fc*#) 

fl,(m, n) = C(m,„)- ^ t? + „ 213+t; 

(10) 


and 


F(m, n) — — — [u(m — 1, n — 1) -j- u(m — 1, n + 1) — ii(m + 1, n — 1) — u(m + 1, n + 1)] 
n 

+ — [u(m — 1, n — 1) — u(m — 1, n + 1) + u(m + 1 , n — 1) — u(m + 1 , n + 1)] 
n 

-J l v ( m ~ I.™ - 1) + v (m - 1 , n + 1) - v(ni + l,n - 1) - v(m + l,n + 1)] 
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1, n — 1) + w{m — 1, n + 1) — w(m + l,n — 1) — w(m + l,n + 1)] (11) 


G(m,n) = — — [6(m — 1, n — 1) + v(m — 1, n + 1) — v(m + l,n — 1) — v(m + 1, n + 1)] 

[#( m ” — 1) - v( m ~ l»n + 1) + 6(m + l,n - 1) - v(m + 1 ,n + 1)] 

[tt(m - 1, n - 1) - u(m - l,n + 1) + u(m + l,n - 1) - u(m + l,n + 1)] 

Jt 

+ -^ [t£(m - l,n - 1) - w(m - l,n + 1) + w(m + l,n - 1) - w(m + 1, n + 1)] (12) 


H(m,n) = — — [t6(m- l,n- l) + w(m-l,n + l)-i6(m + l,n-l)-u>(m + l,n + l)] 

+ ^ [u>(m — l,n — 1) — w(m — l,n + 1) + w(m + l,n — 1) — w(m + l,n + 1)] (13) 
and, everywhere, u>(m,n) is given by 

x m £i(m,n) + n 0(m,n) 

w(m,n) = . (14) 

Since xb can be eliminated from R x and R y by means of the last equation, Eqs. (7) 
and (8) are a system of linear equations for ti(m,n) and 6(m,n) for— oc < m, n < oo. 

The conditions that the bounding box be impenetrable requires that u 9 be zero 
for x = 0 and L\ for all y and z, and v 9 be zero for y = 0 and L 2 for all x and 
z. A consistent set of symmetry conditions which ensure these constraints are the 
following, for all integer values of m,n, 


u(— m, n) = — u(m,n) (15) 

v(—m,n) = v{rri,n) (16) 

u>( — m,n) = ti)(m,n) (17) 

u(m, — n) = 6(m, n) (18) 

v(m, — n) = — 6(m,n) (18) 

w(m, —n) = iy(m,n). (19) 


It follows from these that u(0, n) = 0(m,0) = 0 and that 6(0,0) = 6(0,0) = 0. 
It can be shown that these conditions persist if they are satisfied by the initial 
conditions. These conditions mean that Eqs. (7) and (8) only have to be solved 
on the first quadrant of the m, n plane, i.e., m = 0, 1, 2, . . . , N, n = 0, 1, 2, . . . , N, 
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where N is a large integer where we truncate the system. For instance, if we write 
Eq. (7) with m — 0 and some n, this would require F(0,n), and from Eq. (11) we 
see that we would need £i(-l,n - 1). From Eq. (15), this is equal to -u(l,n - 1), 
which is within the solution set. In general, the column m = -1 and the row n = -1 
are shifted into the solution set. 

From the structure of the equations, one can see that an equation for mode (rn, n) 
is only coupled to four neighboring modes (m — l,n — l),(m — l,n + l), (m + l,n — 1), 
and (m + 1, n + 1). These are the nearest four on diagonal lines through (m, n), and 
these four are each coupled to the nearest four on its diagonals, and so on. This 
means that the modes can be decomposed into two independent “checkerboard” 
subsets, the “even” modes where m + n is an even integer and the “odd” modes 
where m + n is an odd integer. This was pointed out by Bayly (1989). 

In the numerical work, we have taken E = 2 and set the length and time scales by 
talking the maximum vorticity of the basic flow to be unity (A = 1) and L\ = 5.6, 
L 2 = 2.8. The Reynolds number is v~ l in this dimensionless scheme. One could, 
of course, rescale the box so that one of the sides has unit length and rescale the 
Reynolds number and wavenumbers accordingly. The reason for our particular 
choice was determined by the requirements of the nonlinear code which will be 
described in the next section. 

Equations (7) and (8) were solved by a second order Runge-Kutta method over 
a range of values of the parameters R e and k z . Even or odd modes were excited 
by taking initial conditions in which all the mode amplitudes were zero except 
u(2,0) = 1 to generate even modes or u(l,0) = 1 for odd modes. Energy growth 
rates f} were computed from 


„ ,. ln( energy ) 

/? = lira — i — (20) 

t — >oc t } 

where “energy” is the sum of the squares of all the mode amplitudes. It was nec- 
essary to integrate for a long time to approach an asymptote in this formula. Over 
the range computed, we have found that the cutoff N — 40 was adequate. The 
results are presented in two figures. Figure 2 shows the neutral curve (/? = 0) for 
both even and odd modes up to Reynolds numbers of 2000. The band of unstable 
wavenumbers rapidly expands with increasing Reynolds number. Figure 3 shows 
the growth rate (3 versus wavenumber k z at Reynolds number 2000. Computations 
at Reynolds numbers 200, 500, and 1000 are similar, with a pronounced notch in 
the odd mode curve. The growth rate curves computed by Robinson and Saffman 
(1984) for the strained vortex have a similar shape made up of the union of sepa- 
rate growth rate curves for modes with differing internal structure, which look like 
inverted parabolas. The growth rate at their second peak (multiplying their result 
by 2 to get energy growth rate) is about .3 for E - 2, which is comparable to our 
value of about .24. 

A result for unstable plane waves in the unbounded flow with elliptical streamlines 
may be obtained from Landman and Saffman (1987). The result is 
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FIGURE 2. Neutral stability curves, E = 2: even modes; odd modes. 

Heavy dot is at R e number and lowest wave number used in the computation. 



FIGURE 3. Growth rate versus axial wavenumber at R e = 2000, E = 2: even 

modes; odd modes. Vertical dash is at the lowest wave number used in the 

computation. 

where u?o is the uniform vorticity, ctq is the inviscid amplitude growth rate, and 
K 2 = kl + k* + k^. The factor two is inserted to convert to energy growth rate. 
Now ctq/u)q is a function of e/^/[— ( E 2 — 1 )/(E 2 + 1)] where e is strain rate and 7 
is u?o/2. For e/7 < .7 (E < 2.4), it is approximately 

s-KsO- 1 ®*)- 


(22) 
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For an unstable wave with maximum growth rate, the wavevector makes an angle 
of 60° with the rotation axis; therefore, k^ + k* — 3 k 2 z . Using this and v = l/R e 
gives 


JL _ 2 (— _ 1M) 

LJQ y U>Q R e ' 


(23) 


For our case with E = 2 and loq = 1, we get 


/? = .3 



(24) 


For fixed i? e , this gives a curve similar in shape to figure 2 with a maximum of .3 
at k 2 = 0. A neutral curve may be derived from Eq. (24), namely 

k z = .19 R\ /2 . (25) 


We can check a few values. At R e = 2000, 1000, 500, and 200 the formula gives 
8.50, 6.01, 4.25, and 2.68, which compare with the computed values 8.32, 5.81, 3.92, 
and 2.25. The agreement is quite good at the larger Reynolds numbers. Equations 
(24) and (25) should allow a reasonable extrapolation to Reynolds numbers larger 
than those for which we have computed. 

We conclude that the stability results are similar to those for bending waves on a 
concentrated vortex at small wavenumbers and similar to results for an unbounded 
flow with elliptical streamlines at large wavenumber, with comparable growth rates. 


3. Transition to turbulence 

Nonlinear computations have been performed on the Intel Hypercube (i860) using 
a version of Rogallo’s (1981) box code (the “pencil” code) written by Rogallo for 
this parallel processing machine. This code was designed to solve the Navier-Stokes 
equations in a box with periodic boundary conditions. The box is not required 
to have equal sides, but must have volume (27r) 3 . We will take the sides to be 
2L\ x 2X 2 x X3 with the x and y dimensions double the size of the impenetrable 
box. By using proper symmetry in the initial conditions, we can insure that there 
is no flow across the surfaces x = 0, L\ and y — 0,X 2 . The finite dimension in 
the z direction means that X3 is the longest axial wavelength allowed. Defining 
61 = 7t/Xi,6 2 = ^/T 2 as before and 63 = 27r/i3, the volume constraint makes 
b\ 6 2 63 = 1. The parameter 63 is the smallest non- zero axial wavenumber. Allowed 
values are k z = b$k, where k is an integer. Solutions in this periodic system are 
thus of the form 

N 

u{x,y,z,t) = E u(m, n, k) exp i{mb\x + nb^y + kb$z) (26) 

m,n, Jt = — N 


N 


v(x, y, z, t) = Y, i ?(m,n, k)expi(mb\x + nfr 2 y + kb$z) 




(27) 
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FIGURE 4. Surface of 3/8 max vorticity at T = 55.66. Max vorticity is 2.01. 
The tick marks are at .6 intervals in the x-direction and .4 intervals in the y- and 
z-direction. 


N 

w(x,y,z,t ) = E w(rn,n,k)expi(mbix + nb 2 y + kb^z) (28) 

m,n,k= — N 

Conjugate symmetry (tt( — m, — n, — A:) = u(m,n,fc)* , etc.) is imposed to ensure 
that the velocity components are real. Symmetries in the initial conditions like 
those in Eqs. (15)— (19) for each k ensure impenetrability of the side walls. For 
instance, the velocity field of the basic flow, U = 62 sin b\x cos 62 2/ / ( ^ 1 + ^ 2 )y ^ “ 
“61 cos 61 x sin b 2 yj(b\ + 6|) ls generated by the four modes 

6(1, 1,0) = 6(1, -1,0) = —6(— 1,1,0) = — 6(— 1, —1,0) = — 

6(1, 1, 0) = -6(1, -1,0) = — 6(— 1, 1, 0) = 6(— 1, -1,0) = - 4(6 ^ - 6 2 - j (29) 

which have all these symmetries. Perturbations from this can be generated in many 
different ways and will generate different flows. For instance, we could use random 
small initial perturbations. For the flow computed here, we excited both even and 
odd modes by 

u* = .001 sin cos 6362 + .001 sin26jx cos ^Arz 
with v* — 0 and the corresponding w* determined from continuity, that is, 

£t(l, 0, k) = —u{— 1, 0, k) = £(1,0, —k) = —u(— 1, 0, —k) = — .00l2‘/4 


( 30 ) 
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Figure 5. Surface of 3/8 max vorticity at T = 61.71. Max vorticity is 2,78. See 
Figure 4 for the value of the intervals on the axes. 

«(2, 0, k) = -it{~ 2 , 0, k) = fi(2, 0, -*) = — fi(— 2, 0, -Jfc) = -.001i/4 (31) 

for all integer k. With the cosine dependence in the z direction, the endwalls and 
the midplane z = £3/2 will also be impenetrable. 

Finally, we chose to do the computation with E = 2 = L 1 /L 2 and £3 = y/L \ £2 • 
The volume constraint then makes L\ = 5.60, L 2 = 2.8, £3 = 3.96, and then 61 = 
•56, b 2 = 1.12, and 63 = 1 . 59 . We have computed with R t = 2000 and N = 128. 
The lowest wave number in the initial conditions, k z = 1.59, is shown in the stability 
diagrams as a circle in figure 2 and as a vertical dashed line in figure 3. It is near 
the position of maximum growth rate and was chosen for this reason. 

The major results of the computations are presented in figures 4-9, where surfaces 
of constant magnitude of the vorticity axe shown at 6 different times during the 
evolution of this flow. In each figure, the surface is of the vorticity value which is 
3/8 of the maximum vorticity at this time. Both the time and the value of the 
maximum vorticity are given in the figure legends. The coordinate axes shown have 
x to the left, y vertical, and z off to the right. The viewpoint and illumination 
are from the origin. At the initial time (not shown), the surface is approximately 
an elliptic cylinder. The instability causes distortion into a wave oriented roughly 
along a 45° plane with the vortex moving downward and toward positive x at the 
ends of the box and upward and toward negative x in the middle. The secondary 
flow which causes this motion also causes the vortex to be greatly distorted into 
sheet-like structures which resemble the “cups” found by Rogers and Moser (1992) 
in their study of the development of three-dimensional structure in a mixing layer 
after the primary rollup into spanwise rollers. In figures 10a and 10b, we look 
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Figure 6. Surface of 3/8 max vorticity at T = 70.43. Max vorticity is 5.50. See 
Figure 4 for the value of the intervals on the axes. 



FIGURE 7. Surface of 3/8 max vorticity at T = 79.09. Max vorticity is 9.61. See 
Figure 4 for the value of the intervals on the axes. 
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FIGURE 8. Surface of 3/8 max vorticity at T = 83.94. Max vorticity is 13.19. See 
Figure 4 for the value of the intervals on the axes. 



FIGURE 9. Surface of 3/8 max vorticity at T = 87.95. Max vorticity is 15.70. See 
Figure 4 for the value of the intervals on the axes. 
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FIGURE 10. Vorticity contours at the time corresponding to figure 6; a) on the 
plane z = £3/2; b) on the planes 2 = 0 and z = £3. See Figure 1 for the value of 
the intervals on the axes. 



FIGURE 11. Contours of axial velocity on the plane z = £3 /4 at the time 
corresponding to figure 6. Solid contours are positive velocity, dotted contours 
are negative velocity. Maximum velocities are ±.35. See Figure 1 for the value of 
the intervals on the axes. 

at vorticity contours on cuts through the structure shown in figure 6 in order to 
elucidate its structure. In figure 10a, we have taken the cut through the middle of 
the structure along the plane z = £3/2. By symmetry only, the z component of 
vorticity is nonzero on this plane. The maximum vorticity in the box, 5.5, occurs 
at the center of the rather strong round vortex which is apparent here. In the 
three-dimensional rendering this vortex is seen from in front and appears to be like 
a horseshoe vortex. The heavier contour line in figure 10a is at the same vorticity 
level as the surface in the three-dimensional view. In figure 10b, the cut is at the 
ends of the box at planes 2 = 0 and z = £3, where again only ui 2 is non-zero. The 
maximum vorticity in the strong sheet-like structure is 3.5. 

In figure 11, we show contours of the axial velocity w in the plane 2 = £3/4 
through the same structure shown in figure 6. This is more regular than we expected 
and shows characteristics of the 1,1 mode in Robinson and Saffman (1984). The 
first 1 refers to the angular wave number and the second to the number of nodes in 
the axial velocity. Here, despite the already complicated vorticity, the underlying 
secondary flow of the instability is still evident, much amplified, with maximum 
velocities ±.35. Since w is zero at both ends of the box and in the middle, and w on 
the plane 2 = is just the reverse of that shown in figure 11, we can picture the 

secondary flow as two eddies, one in each half of the box, with the velocity toward 
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Figure 12. Energy spectrum at T = 87.95. 

the center of the box in the lower part and away from the center above. Completion 
of this circulation is such as to smash the vorticity downwards and toward positive 
x at the ends while stretching it, and upward and toward negative x in the middle, 
stretching it here also. What is remarkable is that this secondary motion is evident 
throughout the entire sequence of views, even at the time of figure 9 where it is 
quite irregular but with the main trends described above. 

As we proceed from figure 6 to figure 7, the strong vortex which was evident in 
figure 6 has disappeared. The central part has apparently been carried downward 
into the compressive part of the circulation which decreased its vorticity below the 
plotting level. The remnants of its sides have been stretched into the shield-shaped 
vortex sheets on each side of the center. In the next view, figure 8, there is no 
vorticity at the level plotted across the plane z = £ 3 / 2 . This is not clear from 
the figure, but when we rotated it we had an open view down the middle. The 
structures we see are mostly tubelike now, and this is even more evident in figure 
9, which is the time at which the maximum vorticity in the box is largest. At later 
times, this begins to decrease, but the tubelike nature of the vorticity is evident 
still at the time of 103.8 when we ended the computation. 

Figure 12 shows the three dimensional energy spectrum at the time of figure 9. 
While this is quite broad, indicating turbulence, we do not see a —5/3 range at this 
Reynolds number. 

4. Discussion 

Our overall impression of this flow is of two large counter-rotating turbulent 
eddies, each carrying many intense interacting vortex tubes. One can see the de- 
velopment of this pretty clearly in figures 7 and 8, less clearly in figure 9. This 
motion was already present as the growing secondary flow of the instability and 
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FIGURE 13. Surface of 3/8 max filtered vorticity at T = 87.95. Max filtered 
vorticity is 1.19. See Figure 4 for the value of the intervals on the axes. 

was the controlling factor in the development of the vortex sheets. To understand 
the large scale swirling motion without the masking effect of the small vortices, we 
have filtered out the high wavenumber part of the velocity with a Gaussian filter 
with width .2. We have computed the magnitude of the vorticity in this filtered 
flow at the same time as in figure 9. The result is presented in figure 13 as a surface 
of 3/8 of the maximum vorticity, which is 1.19 in this filtered flow. This looks like 
a horseshoe vortex with two vortices wrapping around each of its counterrotating 
legs. The head of the vortex, which is lifted up, contains the point of maximum 
vorticity. The sense of motion is upward between the legs of the vortex and forward 
over the top of the nose, which is the sense of rotation of the original swirling flow. 

The evolution of the sheets into tubes appears to be through the Kelvin- Helmhotz 
instability, or by the Lin-Corcos (1984) version which includes vortex stretching, 
but nothing very dramatic occurs. Local thick spots gradually become thicker and 
develop into tubes. In some preliminary computations at higher Reynolds number, 
the cup- shaped vortex sheet developed into three clearly defined parallel vortex 
tubes. 

In the future, we plan to compute this flow at higher Reynolds number and greater 
resolution with the objective of looking more carefully at the development of small 
scale structure. As a generator of small scale turbulence, this flow is similar to the 
Taylor-Green flow (Brachet et al . , 1983). The present flow has some advantages for 
the study of turbulence because of its development from a typical instability and the 
fact that the turbulent region tends to avoid close contact with the impenetrable 
boundaries. 
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